Phase Transitions in the Distribution of Bipartite Entanglement of a Random Pure State 
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Using a Coulomb gas method, we compute analytically the probability distribution of the Renyi entropies (a 
standard measure of entanglement) for a random pure state of a large bipartite quantum system. We show that, 
for any order q > 1 of the Renyi entropy, there are two critical values at which the entropy's probability distri- 
bution changes shape. These critical points correspond to two different transitions in the corresponding charge 
density of the Coulomb gas : the disappearance of an integrable singularity at the origin and the detachement 
of a single-charge drop from the continuum sea of all the other charges. These transitions respectively con- 
trol the left and right tails of the entropy's probability distribution, as verified also by Monte Carlo numerical 
simulations of the Coulomb gas equilibrium dynamics. 
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Entanglement is a crucial resource in quantum information 
and computation [1] as a measure of nonclassical correla- 
tions between different parts of a quantum system. To ex- 
ploit such correlations to the maximum advantage in quantum 
algorithms, it is desirable to create states with large entan- 
glement. A potential candidate for such a state is a bipartite 
random pure state, its average entanglement entropy being al- 
most maximal [2, 3]. Such a random pure state can also be 
used as a null model or reference point to which the entan- 
glement of an arbitrary time-evolving quantum state may be 
compared. In addition, such random states are also relevant in 
quantum chaotic or non-integrable systems [4, 5]. 

The conclusion that a random pure state in a bipartite sys- 
tem has near maximal entropy is based only on the result for 
the average entropy [2, 3]. Even though the average entropy 
of the random state may be close to its maximal value, the 
probability of the closeness may actually be very small (see 
below). The quantitative evaluation of this probability re- 
quires to compute the full probability distribution of the en- 
tropy, namely its large deviation tails. In addition, the dis- 
tribution of bipartite entanglement may also be used to char- 
acterize entanglement in a multipartite system [6, 7]. In this 
Letter we compute the full distribution of the Renyi entan- 
glement entropies (defined later) for random pure states of a 
bipartite system. The calculation is realized using a Coulomb 
gas method and is valid in the limit when both subsystems are 
large. A by-product of our results is the behavior of the prob- 
ability that the entropy approaches its maximal value In N. 

We start with a standard bipartite system B composed 
of two smaller subsystems A and B, whose respective Hilbert 



spaces H^^^ and Ti.)^"' have dimensions N and M. For sim- 
plicity, we focus on the N = M case, though all our results 
can be extended to the N ^ M case. Let \^) be a nor- 
malized pure state of the full system with its density matrix 
p = IV') ('01 satisfying Tt[p] = 1. The two reduced den- 
sity matrices are denoted pA ~ Trsip] and ps = TiaIp]- 
It is not difficult to prove that both pA and ps share the 
same set of non-negative eigenvalues {Ai, A2, . . . , Ajv} with 
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IZiLi = 1- Let |A^) and |Af ) denote the respective eigen- 
vectors of PA and pb- In this so-called Schmidt basis, an ar- 
bitrary pure state can be represented as 
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lV'>-EVA;|Af)®|A: 



(1) 



This representation is very useful for characterizing the en- 
tanglement between A and B. For example, consider two 
opposite limiting situations: (i) One of the eigenvalues, say 
Ai, is unity and the remaining — 1 are identically zero. 
Then, l^p) = y/Xi (g) |A,f ) factorizes and the system 
is completely unentangled. (ii) All eigenvalues are equal 
(Ai = 1/A^ for all i). Then, all the states are equally present 
in Eq. (1) and the state |0) is maximally entangled. A stan- 
dard measure of entanglement is the von Neumann entropy, 
SvN — — X^i In Ai, which takes its minimum value in 
situation (i) and its maximal value IniV in situation (ii). An- 
other useful measure of entanglement is provided by Renyi's 
entropies, the quantities of major interest here : 
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which also attain their minimum value in situation (i) and 
their maximum value IniV in (ii). As g — 1 and g — > 00, the 
Renyi entropy tends respectively to the von Neumann entropy 
SvN and — In Amax; where A^ax is the largest eigenvalue. 

A pure state ji/;) is called random when it is sampled ac- 
cording to the uniform Haar measure (the unique unitarily 
invariant measure) over the full Hilbert space. As a result, 
the eigenvalues {A^j's also become random variables with the 
joint distribution (for M = N) [3] : 
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Here, (3 — 2 and the (5-function enforces the unit trace con- 
straint Tr[/9A] = 1- Apart from this constraint, (3) is identi- 
cal to the eigenvalue distribution of random Gaussian Wishart 
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(covariance) matrices. For random matrices, the Dyson index 
/3 takes the values 1, 2 or 4 depending on whether the ma- 
trix is real, complex or quaternion. Hence, we shall study (3) 
for general /3, even though /? = 2 in the quantum context. 
The normalization constant Zo can be computed exactly us- 
ing Selberg's integrals [8] as Zq ~ e^^^^ to leading order 
iniV. 

Since the A; 's are random variables distributed as in (3), the 
von Neumann and the Renyi entropies in (2) are also random 
variables. Statistical properties of these observables, as well 
as others such as concurrence, purity, minimum eigenvalue 
etc., have been studied extensively [2, 3, 6, 7, 8, 9, 10, 11]. 
In particular, the average von Neumann entropy (Svn) = 
In TV — 1/2, is close for large N to its maximum value [3]. 
In contrast, few studies have addressed the full distribution of 
the entropy, an exception being the purity E2 = X^il^i '■ 
for small N, the distribution of purity is known exactly ['^J; 
for large N, the Laplace transform of the purity distribution 
was studied recently for positive Laplace variables [7]. How- 
ever, the inverse Laplace transform of this quantity provides 
only partial information about the purity distribution. 

The goal of our Letter is to compute analytically, for large 
N and all g > 1, the full distribution of the Renyi entropies 
in (2), or equivalently of ~ J2iLi K ~ cxp[(l — q)Sq]. 
The quantities Sg satisfy the inequalities N^""^ < Eg < 1 for 
q > 1, with the upper and lower bounds corresponding to the 
unentangled (i) and the maximally entangled (ii) situations. 
The distribution of E, is written using (3) as 



P(E, = W'-' s) 



(4) 



The approach we employed to treat (4) is a saddle-point 
method to identify the configuration of the eigenvalues {A^j's 
that dominates for large A^. Configurations at large TV 
are characterized by the continuous density p(A, N) = 
TV~^^ji5(A — Xi) and the main challenge, accomplished 
here, is to find the saddle-point density pdX, N). 

Let us first summarize our main results. The normalization 
^iLi ~ 1 implies that the typical amplitude of the eigen- 
values Xi ^ 1 /N and hence E^ ~ N^^i for large N. We de- 
fine then the rescaled intensive variable s = N^^^ E^ whose 
lowest value s = 1 corresponds to the maximally entangled 
situation (ii). In Fig. 1, a typical plot of P(S, = TV^'^s, TV) 
vs s is shown for large TV and fixed q > 1: the distribution 
has a Gaussian peak flanked on both sides by non-Gaussian 
tails. Specifically, we find two critical values s = si{q) and 
s ~ s{q) separating three regimes I (1 < s < si ((?)), II 
(siiq) < s < s{q)) and III (s > s{q)). At the first criti- 
cal point si{q) the distribution has a weak singularity (third 
derivative is discontinuous). At the second critical point s{q), 
a Bose-Einstein type condensation transition occurs (see be- 
low). These changes are a direct consequence of two phase 
transitions in the associated optimal charge density (shown in 
Fig. 2). In regime I, the optimal charge density has a compact 
support [(^i/N,C2/N], where ^1 is strictly positive (Fig. 2a). 




FIG. 1 : The schematic distribution of E,, 
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as a function of s for fixed large TV. Two critical points s = si{q) 
and s = s{q) separate three regimes I, II and /// characterized 
by the different optimal densities shown in Fig. 2. The maximally 
entangled state s = 1 is at the extreme left, in the large deviation 
tail well-spaced from the average s{q). (Inset) The large deviation 
functions $ for the distribution of Eq, in the three different regimes. 
Analytical predictions (red solid line) are compared to the results 
(blue points) of Monte Carlo numerical simulations of the Coulomb 
gas equilibrium dynamics. 
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FIG. 2: Scheme of the optimal saddle-point density pc of the eigen- 
values, or equivalently of the Coulomb gas charges, for 1 < s < 
si(q) (regime I), si(g) < s < s{q) (II) and s > s{q) (III). In the 
regime III, the maximum eigenvalue becomes larger than all the oth- 
ers (~ 0(1/TV)), as shown by the isolated bump at t in (2c). 



When the control parameter s exceeds si{q) (regime II), the 
left edge of the support sticks to zero while the upper edge 
L /TV moves to the right as s increases (Fig. 2b), till the sec- 
ond critical value s{q) corresponding to L = 4. For s > s{q) 
(regime III), we find that one eigenvalue (the small bump in 
Fig. 2c) splits off the sea of all the other TV — 1 eigenvalues, 
which remain ^ 0(1/TV). This second phase transition is 
reminiscent of the real-space condensation phenomenon ob- 
served in a class of lattice models for mass transport, where a 
single lattice site carries a thermodynamically large mass [12]. 

Note that for q = 2, the presence of two phase transitions 
was also noticed in Ref. [7] for the Laplace transform of the 
distribution P(E2,TV). However, the nature of the regime 
III was not elucidated and the corresponding optimal den- 
sity and the partition function were not calculated. To derive 
the full distribution P(E2, TV), one needs the partition func- 
tion in all three regimes, which is what we do here for all 
q > 1 at large TV. We also find exact expressions for the two 
critical points: A-is^q) = r(g + l/2)/(77iT(g + 2)) and 
(4(<z + l)/3g)-9si(g) = r(g + 3/2)/(0Fr(g + 2)). From 
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the Gaussian form near the peak s = s{q), we also read off 
the mean and the variance of the entropy Sq for all q : 



(5,) «ln(7V)-i^; Var(5,) 



(5) 



Let us now briefly outline how the previous results are de- 
rived. By using (3) and (4), we obtain 

P(I]„ N) = ^ ; Z(E,) = I e-^^«^-» d\. , (6) 

with E{{X,}) = -(1/2- ln(AO - E,<, 1^ |A. - I 
and the integral runs over the subspace satisfying the two con- 
straints, J2i Ai = 1 and J2i K — ^9- The expression for 
E{{Xi}) is interpreted as the energy of a Coulomb gas of 
charged particles with coordinates Ai that repel each other via 
2-d logarithmic interactions and are also subject to an exter- 
nal logarithmic potential. In the large limit, we can char- 
acterize the configuration of the Coulomb gas' charges by the 
normalized density p{\, N) = S{\ — A;). Due to the 



constraint Ei A; = 1, typically A; ^ 1/N. Hence, the charge 
density scales as p(A, N) w N p{XN) and we introduce the 
rescaled variable s = N''^^T^q. We then replace the mul- 
tiple integral in (4) by a functional integral over all possible 
normalized and rescaled charge density functions p{x) satis- 
fying the three constraints; J p{x)dx = 1, J xp{x)dx = 1 
and J x'^p{x)dx = s. The resulting functional integral over 
p{x) is evaluated in the large N limit via the saddle point 
method. This constrained Coulomb gas approach has been 
used successfully in a variety of contexts that include the dis- 
tribution of the top eigenvalues of Gaussian and Wishart ma- 
trices [13, 14, 15], phase transition in the restricted trace en- 
semble [ 1 6], purity partition function in bipartite systems [7], 
nonintersecting Brownian interfaces [ ], quantum transport 
in chaotic cavities [18], information and communication sys- 
tems [19], and the index distribution for Gaussian random 
fields [20, 2 1 ] and Gaussian matrices [22]. 



The constrained Coulomb gas approach yields 

oc e-'S^'^^M. To the lead- 



ing order in N, the effective energy reads 
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dx dx' p{x)p{x')\n\x ~ x'\ +/^o( / dxp{x) — 1) +pi[ / dxxp{x) — 1) + / dxx'^ p{x) — s) ,(7) 



where the Lagrange multipliers pQ, pi and p2 enforce the con- 
straints. For large N, the method of steepest descent gives: 

P(l]g = iV^"«s) cx £"'^^'^=['''^1 where pc{x) minimizes 
the energy : SEg [p]/5p = 0. This gives the integral equation 

/>oo 

V{x) = pq + pix + p2x'^ ~ / pc{x') In \x — x'\dx' , (8) 

with V{x) acting like an effective external potential. Differ- 
entiating once more with respect to x leads to 

p,+qp^x^-^ =V r ^^dx', (9) 

Jq X X 

where V denotes Cauchy's principal part. The single-support 
solution to (9) is found by using Tricomi formula [23] and 
yields the regimes sketched in Figs. 1 and 2. 

Regime I: For 1 < s < si(g), we find that pi < 0, p2 > 
and the effective potential V{x) has a minimum at a nonzero 
X. This indicates that the charges concentrate around this 
nonzero minimum over a support [Ci, (^2] for all g > 1 (see 
Fig. 2a). For q ~ 2, the edges C1.2 = 1 =F 2-y/s — 1 and 
the solution pc{x) = \J {(2 — x){x — (i)/{2tt{s — 1)) van- 
ishes at both edges. This solution exists for Ci > 0, i.e., for 
s < si{2) = 5/4, and the distiibution P{Y.2 = s/N,N) ^ 
g-PN <i>i(s) ^jjjj jjjg j^gg deviation function 

ci>,(s) = -iln(s-l)-i. (10) 



The behavior for g 7^ 2 is qualitatively similar, though the ex- 
pressions are cumbersome [24]. Setting s = 1 + e around the 
maximal entropy state s = 1, the probability at this extreme 
left tail scales as ~ i.e. it is very small for large N. As 

s approaches si{q) from below, pi and the minimum of V{x) 
tend to zero. This signals that the charges now concentrate 
near the origin and the onset of regime II. 

Regime II: For si{q) < s < s{q), the charges concentrate 
over a support [0, L] (see Fig. 2b). For q ~ 2, the optimal 
charge density takes the simple form Pc{x) ~ ^ (L — a;)(A + 
Bx)/{Tr^), where A = 4(L - 2)/L ^ S = 8(4 - L)/L^ 
and the right edge L = 6 — 2\/9 — 4s. Evaluating the energy 
for large N, we get ^(Sa = s/N, N) - e-'3^'*"(-^) with 

*//(«) = -^ln(i/4) + ;^-^ + ^. (11) 

Comparing (10) and (1 1), it is verified that the large deviation 
functions match at the critical point si — 5/4 up to the second 
derivative, while the third is discontinuous: $} (5/4) = —32 
and $jj''(5/4) = —16. The function <I'//(s) is quadratic 
$77(5) « (s — 2)^/8 around its minimum at s = s(2) = 2. 
Thus, the distribution P(I]2 = s/N, N) has a Gaussian peak 
near s = 2, with the mean {Y,2) = 2/N and the variance 
Var(S2) = 4/ {/3N'^) for large N. The corresponding expres- 
sions for arbitrary q > 1 are given in (5). 

For any q > 1, pi is positive, /i2 — > as s s{q) and 
P2 < for s > s{q). This indicates that the potential V{x) in 
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(8) becomes non-monotonic for s > s{q) : it increases around 
the origin, reaches a maximum at x* — {—^i/ qpL2)'^^^'^~^^ 
and then decreases monotonically for x > x* . It follows that 
the minimum at a; = as well as the associated saddle-point 
solution become metastable. This solution still exists over a 
finite range for s > s{q) (e.g., for q = 2 over 2 < s < 9/4). 
For s > s{q), however, there is a second stable solution where 
one eigenvalue splits off the sea of the remaining {N — 1) 
eigenvalues (see regime III below). For s > s{q), the energy 
associated with this stable solution is lower by ~ 0{1/N) 
only, as compared to the energy of the metastable solution. 

Regime III: For s > s{q), the correct density of states 
consists of two disjoint parts: (a) {N — 1) eigenvalues re- 
main ^ 0{1/N), concentrated over a finite support includ- 
ing the origin; (b) the top eigenvalue Amax takes a larger 
value and moves away from the sea of all the others (see 
Fig. 2c). The saddle point method thus needs to be slightly 
revised. For simplicity, let us focus only on g = 2. We write 
Amax = t and label the remaining (TV — 1) eigenvalues by 
their continuous density p(A) = J2i^max ^i)- We 
then express the energy ^^[{Ai}] in (6) in terms of p{X) and 
t, treating both of them as variables. This gives P{T^2 = 
S,N) ex / Vp J dt e~^-^^^P'*\ where the effective energy 
Hs [p, t] has a long expression that includes p, t and three La- 
grange multipliers enforcing the constraints [ ]. Assuming 
that p(A) has a finite support over [0, Q with C < i, we min- 
imize the effective energy over both p and t. The equations 
SHs/6p = = dHs/dt are solved again using Tricomi's 
theorem [23]. Substituting the solutions for p{X) and t in 
the effective energy finally yields the distribution P(S2, N) 
at the leading order in N. We have verified that in the regime 
2 < s < 9/4, the resulting distribution coincides with that of 
regime II, i.e. the transition at E2 2/N is smooth. The 
maximum eigenvalue Amax = t dominates at the upper edge 
of the regime III, when S2 ~ 0(1), and we find [24] that 

P(S2 = ^, iV) ~ (1 - 

Numerical Simulations. To verify analytical predictions, 
we simulated the distribution (3) of the eigenvalues A;, which 
is interpreted as the Boltzmann weight of a Coulomb gas 
and sampled using the Metropolis algorithm (see, e.g., [25]). 
Specifically, we start with a configuration of the A^'s satisfy- 
ing J2i ~ 1- The moves in the Metropolis scheme consist 
of picking at random a pair (Ai, Xj ) and proposing to modify 
them as (A^ + e, Xj — e) where e is set to achieve the standard 
average rejection rate 1/2 [25] . As usual, the move is accepted 
with probability e"'^^^ if AE > and with probability 1 if 
AE < 0, where AE is the change in energy i?[{Ai}] (the 
move is rejected if one of the eigenvalues becomes negative). 
This ensures that at long times we reach thermal equilibrium 
with the correct Boltzmann weight cx e^^-^^^^^^^ satisfying 
the constraint = 1- We then construct the histogram of 

P(Eq — J2i Numerical data compare very well with 

our analytical predictions (see the inset of Fig. 1 for q — 2, 
N ~ 50) ; we also verified that a single eigenvalue detaches 
from the sea in regime III (intuitively, multiple drops are unfa- 
vorable as they compress the sea more than a single drop due 



to the convexity of ^ A' for q > 1). 

In conclusion, we have obtained the first complete charac- 
terization of the quantum entanglement's statistical properties 
in a bipartite random pure state of large dimensions N. The 
average of the Renyi entropies is indeed close to its largest 
value IniV. This is, however, the mere consequence of the 
typical amplitude 1/iV of the density matrix eigenvalues. The 
distribution of the eigenvalues mostly affects the 0(1) con- 
tribution to the entropy. The probability to approach In N is 
actually found to decay rapidly at large N, as clearly shown by 
the full probability distribution of the Renyi entropies derived 
here. The spreading of the eigenvalues becomes prominent 
in the regime III (in Figs. 1 and 2) where a condensation oc- 
curs and the contribution by the single top eigenvalue of the 
density matrix is thermodynamically relevant. 
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